home *** CD-ROM | disk | FTP | other *** search
/ C/C++ Users Group Library 1996 July / C-C++ Users Group Library July 1996.iso / vol_100 / 157_01 / fft.c < prev    next >
Text File  |  1987-10-10  |  5KB  |  205 lines
  1. /* 
  2. HEADER:     CUG
  3. TITLE:        FFT.C - Fast Fourier Transform
  4. VERSION:    1.00
  5. DATE:        05/27/85
  6. DESCRIPTION:    "A Fast Fourier Transform implementation based on
  7.         Cooley's successive-doubling method."
  8. KEYWORDS:    fft, filter, fourier, transform
  9. SYSTEM:        Any
  10. FILENAME:    FFT.C
  11. WARNINGS:    "Complex numbers are represented by their real
  12.         and imaginary components in a 2-D array. Data
  13.         must be presented in multiples of two."
  14. CRC:        xxxx
  15. SEE-ALSO:    FWT.C
  16. AUTHORS:    Ian Ashdown - byHeart Software
  17. COMPILERS:    Any C compiler
  18. REFERENCES:    AUTHORS: R.F. Gonzalez & P. Wintz;
  19.         TITLE:     Digital Image Processing;
  20.         AUTHORS: J. Cooley, P. Lewis & P. Welch;
  21.         TITLE:     The Fast Fourier Transform and Its
  22.              Applications
  23.              IEEE Trans. Educ.
  24.              pp. 27 - 34, Vol. E-12, No. 1, 1969;
  25.         AUTHORS: A. Aho, J. Hopcroft & J. Ullman;
  26.         TITLE:     The Design and Analysis of Computer
  27.              Algorithms
  28.              Chapter 7, The Fast Fourier Transform
  29.              and Its Applications;
  30. ENDREF
  31. */
  32.  
  33. /*-------------------------------------------------------------*/
  34.  
  35. #include "math.h"
  36.  
  37. typedef struct complex    /* Complex number structure */
  38. {
  39.   double real,
  40.      imag;
  41. } COMPLEX;
  42.  
  43. #define PI  3.1415926536
  44.  
  45. /* FFT() - Fast Fourier Transform accepts a two dimensional
  46.  *       "double" array of two rows of N+1 elements, where the
  47.  *       zeroth column elements are not used and N is equal to
  48.  *       an integer power of two, and a variable n, set equal
  49.  *       to N above. The discrete Fourier transform is computed
  50.  *       in place and returned in the same array.
  51.  */
  52.  
  53. void fft(n,point)
  54. int n;
  55. COMPLEX point[];
  56. {
  57.   register int a,
  58.            b;
  59.   int d = 1,
  60.       e,
  61.       f;
  62.   COMPLEX u,
  63.       w,
  64.       temp;
  65.   void c_bitrev();
  66.  
  67.   /* Reorder complex data vector by bit reversal rule */
  68.  
  69.   c_bitrev(n,point);
  70.  
  71.   /* Perform successive doubling calculations */
  72.  
  73.   while(d < n)
  74.   {
  75.     e = d;
  76.     d *= 2;
  77.     u.real = 1.0;
  78.     u.imag = 0.0;
  79.     w.real = cos(PI/e);
  80.     w.imag = -sin(PI/e);
  81.     for(b = 1; b <= e; b++)
  82.     {
  83.       for(a = b; a <= n; a += d)
  84.       {
  85.     f = a + e;
  86.     temp.real = point[f].real * u.real -
  87.         point[f].imag * u.imag;
  88.     temp.imag = point[f].real * u.imag +
  89.         point[f].imag * u.real;
  90.     point[f].real = point[a].real - temp.real;
  91.     point[f].imag = point[a].imag - temp.imag;
  92.     point[a].real = point[a].real + temp.real;
  93.     point[a].imag = point[a].imag + temp.imag;
  94.       }
  95.       temp.real = u.real;
  96.       u.real = u.real * w.real - u.imag * w.imag;
  97.       u.imag = temp.real * w.imag + u.imag * w.real;
  98.     }
  99.   }
  100.  
  101.   /* Normalize transformed data vector */
  102.  
  103.   for(a = 1; a <= n; a++)
  104.   {
  105.     point[a].real = point[a].real/n;
  106.     point[a].imag = point[a].imag/n;
  107.   }
  108. }
  109.  
  110. /* IFFT() - Inverse Fast Fourier Transform accepts a two
  111.  *        dimensional "double" array of two rows of N+1
  112.  *        elements, where the zeroth column elements are not
  113.  *        used and N is equal to an integer power of two, and a
  114.  *        variable n, set equal to N above. The inverse
  115.  *        discrete Fourier transform is computed in place and
  116.  *        returned in the same array. 
  117.  */
  118.  
  119. void ifft(n,point)
  120. int n;
  121. COMPLEX point[];
  122. {
  123.   register int a,
  124.            b;
  125.   int d = 1,
  126.       e,
  127.       f;
  128.   COMPLEX u,
  129.       w,
  130.       temp;
  131.   void c_bitrev();
  132.  
  133.   /* Reorder complex data vector by bit reversal rule */
  134.  
  135.   c_bitrev(n,point);
  136.  
  137.   /* Perform successive doubling calculations */
  138.  
  139.   while(d < n)
  140.   {
  141.     e = d;
  142.     d *= 2;
  143.     u.real = 1.0;
  144.     u.imag = 0.0;
  145.     w.real = cos(PI/e);
  146.     w.imag = sin(PI/e);
  147.     for(b = 1; b <= e; b++)
  148.     {
  149.       for(a = b; a <= n ; a += d)
  150.       {
  151.     f = a + e;
  152.     temp.real = point[f].real * u.real -
  153.         point[f].imag * u.imag;
  154.     temp.imag = point[f].real * u.imag +
  155.         point[f].imag * u.real;
  156.     point[f].real = point[a].real - temp.real;
  157.     point[f].imag = point[a].imag - temp.imag;
  158.     point[a].real = point[a].real + temp.real;
  159.     point[a].imag = point[a].imag + temp.imag;
  160.       }
  161.       temp.real = u.real;
  162.       u.real = u.real * w.real - u.imag * w.imag;
  163.       u.imag = temp.real * w.imag + u.imag * w.real;
  164.     }
  165.   }
  166. }
  167.  
  168. /* C_BITREV() - Reorder complex data vector by bit reversal
  169.  *        rule.
  170.  */ 
  171.  
  172. void c_bitrev(n,point)
  173. int n;
  174. COMPLEX point[];
  175. {
  176.   register int i,
  177.            j,
  178.            k;
  179.   COMPLEX temp;
  180.  
  181.   j = n/2 + 1;
  182.   for(i = 2; i < n; i++)
  183.   {
  184.     if(i < j)
  185.     {
  186.       temp.real = point[i].real;
  187.       temp.imag = point[i].imag;
  188.       point[i].real = point[j].real;
  189.       point[i].imag = point[j].imag;
  190.       point[j].real = temp.real;
  191.       point[j].imag = temp.imag;
  192.     }
  193.     k = n/2;
  194.     while(k < j)
  195.     {
  196.       j = j - k;
  197.       k /= 2;
  198.     }
  199.     j = j + k;
  200.   }
  201. }
  202.  
  203. /* End of FFT.C */
  204. .imag;
  205.     temp.imag = po